Thunderstrom Analysis

Imports

library(evd); 
library(evdbayes);
library(coda);
library(ismev);

Lade nötiges Paket: mgcv
Lade nötiges Paket: nlme
This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.

library(xts);

Lade nötiges Paket: zoo

Attache Paket: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading the Data

load("../data/CAPE_Minder_Rychener_Malsot.RData")
load("../data/NINO34.RData")
load("../data/SRH_Minder_Rychener_Malsot.RData")

Generate PROD

prod = sqrt(cape)*srh

** Create Time Series Objects **

dates <- seq.Date(as.Date("1979-1-1"), as.Date("2015-12-31"), by=1)
feb29ix <- format(as.Date(dates), "%m-%d") == "02-29"
dates <- dates[!feb29ix]
prod_ts = xts(prod, order.by = rep(dates, each=8))

Beginning the Analysis

month_names = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
get_monthly = function(x) {
  output = list()
  len = nrow(x)
  
  # Get month of element
  month = time(x)
  month = gsub("....-", "", month)
  month = gsub("-..", "", month)
  monthlist = unique(month)
  for (i in 1:12) {
    output[[i]] = x[month == monthlist[i],]
  }
  names(output) = month_names
  return(output)
}
monthly_max = get_monthly(apply.monthly(prod_ts, max))
r = 2
r_monthly_max = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:r]))

# get the monthly maxima of prod
m1 = as.data.frame(apply.monthly(prod_ts, max))$V1;
# produce the gumbel qq plot
gumbel_qq = function(x, title="") qqplot(qgumbel(c(1:length(x))/(length(x)+1)),
                                         x,
                                         main = paste("Gumbell Q-Q Plot", title),
                                         xlab = "Theoretical Quantiles", 
                                         ylab = "Sample Quantiles") ;
gumbel_qq(m1)

#qqplot(qgumbel(c(1:length(m1))/(length(m1)+1)),m1,main = "Gumbell Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles") ;

The qq plot doesn’t fit very well, especially in the lower tail. This is likely due to seasonal dependence.

Fitting GEV to the entire data

# fit gevd with MLE and produce diagnostic plots
fitmax.MLE<-fgev(m1,method="Nelder-Mead")
par(mfrow=c(2,2))
fitmax.MLE

Call: fgev(x = m1, method = "Nelder-Mead") 
Deviance: 9272.599 

Estimates
      loc      scale      shape  
7.519e+03  7.013e+03  4.757e-03  

Standard Errors
      loc      scale      shape  
421.90201  331.43749    0.05347  

Optimization Information
  Convergence: successful 
  Function Evaluations: 171 

plot(fitmax.MLE)

Poor fit, probably because the distribution isn’t stationary. This is especially visible in the Probability plot, in which the confidence band is surpassed, indicating a poor fit.

# fit gevd with Bayesian Techniques
# use the previous outputs (rounded) as initial values (use a different shape)
init<-c(1.6e4,4e3,0.1)
# arbitrary priors
mat <- diag(c(10000,10000,100)) 
pn <- prior.norm(mean=c(0,0,0),cov=mat)
# proposal standard deviation (found by trial and error to get 20%<acceptance rate<40%)
psd<-c(500,0.03,0.02)
# draw 3k samples from markov chain
nit = 30000
thinning = 300
post<-posterior(nit, init=init,prior=pn,lh="gev",data=m1,psd=psd)
# diagnostic plots
MCMC<-mcmc(post[1:nit %% thinning == 0, ], thin=300) 
plot(MCMC) 

attr(mcmc(post),"ar")

            mu sigma   xi total
acc.rates 0.24  0.71 0.70  0.55
ext.rates 0.00  0.00 0.02  0.01

#MCMC_stab <- mcmc(post, thin=50, start=1000)
acf(mcmc(post[1:nit %% thinning == 0, ]))

We observe that there seem to be no substantial issues from the autocorrelation plots.

apply(mcmc(post[1:nit %% thinning == 0, ]),2,mean)

           mu         sigma            xi 
  214.8136408 11344.8681059    -0.1940404 

apply(mcmc(post[1:nit %% thinning == 0, ]),2,sd)

          mu        sigma           xi 
100.50420904 741.34912578   0.03950923 

** Fit with MLE for months separately**

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
monthly_fits = list()
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  
  if (i %in% error_cases) {
    monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  }
  else {
    monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead")
  }
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9

optimization may not have succeeded

[1] 10
[1] 11
[1] 12

names(monthly_fits) = names(monthly_max)

gumbel_qq(data.frame(monthly_max[[9]])[,1], "September")

gumbel_qq(data.frame(monthly_max[[12]])[,1], "December")

get_se = function(x, ix) {
  if (is.null(x$std.err)) 0
  else x$std.err[ix]
}
mle_loc = unlist(lapply(monthly_fits, function(x) x$estimate[1]))
mle_loc_se = unlist(lapply(monthly_fits, get_se, 1))
mle_scale = unlist(lapply(monthly_fits, function(x) x$estimate[2]))
mle_scale_se = unlist(lapply(monthly_fits, get_se, 2))
mle_shape = unlist(lapply(monthly_fits, function(x) x$estimate[3]))
mle_shape_se = unlist(lapply(monthly_fits, get_se, 3))

plot_w_err = function(x, y, se, title = NULL) {
  max_ix = which.max(y)
  min_ix = which.min(y)
  plot(x, y,
       ylim = c(y[min_ix] - se[min_ix], y[max_ix] + se[max_ix]),
       main = title)
  arrows(x,y-se,x,y+se, code=3, length=0.02, angle = 90)
}
plot_w_err(1:12, mle_loc, mle_loc_se, "Location Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_scale, mle_scale_se, "Scale Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_shape, mle_shape_se, "Shape Parameter as Estimated with Likelihood")

** Fit with Bayesian Methods for months separately**

# Fits GEV distribution with bayesian method for given parameters
bayes_fitter = function(x, 
                        init = c(1e3, 1e3, 0.1), # Initial values
                        mat = diag(c(10000,10000,100)),
                        psd = c(500,0.1,0.1), # Proposed SDev
                        nit = 3000, # Nb Iterations
                        thinning = 50, # Thinning Factor
                        do_diagn = FALSE, # Bool whether to show diagnostic plots
                        do_autoreg = FALSE, # Bool whether to show autoreg plots
                        seed = 42 # Seed 
                        ) {
  set.seed(seed)                
  pn = prior.norm(mean=c(0,0,0),cov=mat)
  post = posterior(nit, init=init, prior=pn, lh="gev", data=x, psd=psd)
  
  if(do_diagn) {
    MCMC<-mcmc(post[1:nit %% thinning == 0, ]) 
    plot(MCMC) 
  }
  if(do_autoreg) {
    acf(mcmc(post[1:nit %% thinning == 0, ]))
  }
  list(posterior = post, 
       acceptance_rate = attr(mcmc(post),"ar"))
}
# Iteratively fits GEV with bayesian methods, until the fit has 
# acceptable acceptance rates (i.e. 0.2 < AR < 0.4). If the AR is too high, 
# the proposed SDev for the parameter is multiplied with 1.5. If it's too small,
# the proposed SDev is divided by 2. This is repeated until the acceptance rate
# is good for all parameters, or max_it is reached. Then, a final model is fitted
# with more iterations. 
bayes_fitter_param_search = function(x,
                                     psd_init = c(500,0.1,0.1), # Initial proposed SDev
                                     nit_full = 3000, # Nb Iterations for final model
                                     nit_search = 150, # Nb Iterations for param search
                                     do_diagn = FALSE, # Bool whether to show diagnostic plots
                                     do_autoreg = FALSE, # Bool whether to show autoreg plots
                                     max_it = 20,
                                     ... # Additional params passed to bayes_fitter
                                     ) 
{
  # Iterate until desired acceptance rate is obtained
  cont = TRUE
  psd = psd_init
  it = 0
  while(cont) {
    it = it+1
    if (it > max_it) {
      warning("The ")
    }
    fit = bayes_fitter(x, psd=psd, nit=nit_search, do_diagn=FALSE, 
                       do_autoreg=FALSE,...)
    acc_rates = fit$acceptance_rate[1, 1:3]
    
    too_high = acc_rates > .4
    too_low = acc_rates < .2
    
    if (all(!too_high) && all(!too_low)) { # All acceptance rates lie within threshold
      cont=FALSE
    } else if (it > max_it) { # max_it is reached
      cont=FALSE
      warning("max_it was reached")
    } else { # Correct values which have wrong threshold
      psd[too_high] = psd[too_high] * 1.5
      psd[too_low] = psd[too_low] / 2
    }
  }
  
  # Fit final model
  bayes_fitter(x, psd=psd, nit=nit_full, do_diagn=do_diagn, 
               do_autoreg=do_autoreg, ...)
  
}
                                     
monthly_bayes_fit = lapply(monthly_max, bayes_fitter_param_search, do_diagn = TRUE, 
                           do_autoreg = FALSE,
                           psd = c(500,0.3,0.3), nit_full=30000, nit_search = 30000,
                           thinning = 300)

TODO-> R largest fit

PART 2 First, we check if the location parameter depends on time using a likelyhood ratio test

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  
  fit_const = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             nsloc = trend,
                             std.err = FALSE)
  
  ratios[[i]] = fit_const$dev-fit_dependant$dev 
}

[1] 1

optimization may not have succeeded

[1] 2

optimization may not have succeeded

[1] 3
[1] 4

optimization may not have succeeded

[1] 5

optimization may not have succeeded

[1] 6
[1] 7
[1] 8
[1] 9

optimization may not have succeeded

[1] 10
[1] 11

optimization may not have succeeded

[1] 12

names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for time independance, Bonferroni multiple Testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Now, let’s check for independance from ENSO

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
# split nino data into months
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  trend = n[i,]
  trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
  fit_const = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             nsloc = trend,
                             std.err = FALSE)
  
  ratios[[i]] = fit_const$dev-fit_dependant$dev 
}

[1] 1

optimization may not have succeeded

[1] 2

optimization may not have succeeded

[1] 3
[1] 4

optimization may not have succeeded

[1] 5
[1] 6

optimization may not have succeeded

[1] 7
[1] 8

optimization may not have succeeded

[1] 9

optimization may not have succeeded

[1] 10
[1] 11
[1] 12

names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for independance from ENSO, Bonferroni multiple testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Another method is the chi plot:

# plot the chi plot for dependance analysis
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
for (i in 1:length(monthly_max)) {
  nino = n[i,]
  m_data=as.data.frame(monthly_max[[i]])$V1
  dat.m1_month = cbind(m_data,trend);
  dat.m1_nino = cbind(m_data,nino);
  par(mfrow=c(2,2))
  chiplot(dat.m1_month,main1 = "Chi Plot Time",main2 = "Chi Bar Plot Time");
  chiplot(dat.m1_nino,main1 = "Chi Plot ENSO",main2 = "Chi Bar Plot ENSO");
}

PART 3 We will now analyse temporal clustering of extremes. For this, we will use the exiplot function from the evd library.

# define a function for getting the extremal indices for each month for a given threshold
monthly_eindex <- function(data, threshold_p, r=0){
  ei = list()
  for (i in 1:length(data)) {
    threshold = quantile(as.data.frame(data[[i]])$V1, threshold_p)
    ei[[i]]=exi(as.data.frame(data[[i]])$V1, threshold, r)
  }
  names(ei) = names(data)
  return(ei)
}
ei = monthly_eindex(get_monthly(prod_ts), 0.95)
plot(unlist(ei), main="Extremal Index by Month, 95%-Quantile as Threshold", xlab="Month", ylab="Extremal index")

We observe that the extremal index is _0.25-0.45, we can therefore conclude that we have strong dependance of extremes, with the limiting mean cluster size being roughly from 2 to 4. The clustering has no effect for estimaters based on the (monthly) maximum, but the r largest estimator is influenced by it.

PART 4 First, let’s estimate the 10 year return level using point process approach

monthly_fits_pp = list()
monthly_data = get_monthly(prod_ts)
#error_cases = c(1,2,3,4,5,6,7,8,9,10,11,12)
month_days = c(31,28,31,30,31,30,31,31,30,31,30,31)
for (i in 1:length(monthly_max)) {
  print(i)
  threshold = quantile(as.data.frame(monthly_data[[i]])$V1, 0.95)
  
  if (i %in% error_cases) {
    monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
                             threshold = threshold,
                             model="pp",
                             npp = month_days[i]*8,
                             cmax = TRUE,
                             r = 1,
                             std.err = FALSE,
                             method = "Nelder-Mead")
  }
  else {
    monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
                             threshold = threshold,
                             model="pp",
                             npp = month_days[i]*8,
                             cmax = TRUE,
                             r = 1,
                             method = "Nelder-Mead")
  }
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12

names(monthly_fits_pp) = names(monthly_data)
for(i in 1:12){
  par(mfrow=c(2,2)) 
  plot(monthly_fits_pp[[i]])
}

The fit in february has completely failed, and the others are not very good either

We will still estimate the return levels:

return_level = function(x,period=20){
  p = 1/period
  loc = x$estimate[[1]]
  scale = x$estimate[[2]]
  shape = x$estimate[[3]]
  level = loc + scale*(((-log(1-p))^-shape-1)/shape)
  return(level)
}
return_level_20 = lapply(monthly_fits, return_level) # 20 for testing
return_level_100 = lapply(monthly_fits, return_level, period=100)
return_level_50 = lapply(monthly_fits, return_level, period=50)
plot(unlist(return_level_100),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

plot(unlist(return_level_50),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

TODO: estimat with mcmc Assuming that we have the posterior densities of the markov chains, call theis function to plot return level histograms

return_level_mcmc = function(posterior,period=20){
  u = mc.quant(posterior,p=1-1/period,lh="gev")
  label_mcmc_rl = sprintf("%s Year return level",period)
  hist(u,nclass=20,prob=T,xlab=label_mcmc_rl)
}